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A generalization of the lattice Bhatnagar-Gross-Krook (LBGK) model for the simulation of hy- 
drodynamics is presented, which takes into account the difference and the frame-independence of 
^^ ' the relaxation of non-hydrodynamic modes. The present model retains the computationally efficient 

^^ \ standard LBGK form with the generalized equilibrium explicitly derived. The two-dimensional re- 

alization on the standard lattice is discussed in detail. Performance of the model is assessed through 
a shear layer simulation and enhanced stability and accuracy with respect to the standard LBGK 
are reported. The results demonstrate that the present model is a useful upgrade of the standard 
C^^ , LBGK without compromising its computational efficiency and accuracy. 
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INTRODUCTION 



The lattice Bhatnagar-Gross-Krook (LBGK) model was conceived about twenty years ago [l|, Q and has rapidly 
. ' taken a dominant role in the lattice Boltzmann approach to the simulation of complex hydrodynamic phenomena 



@j 13 • The success of LBGK is primarily based on its computational efficiency, accuracy and stability at moderate 

Reynolds number simulations. However the standard LBGK has its limits in addressing direct numerical simulation of 

high Reynolds number flows caused by severe numerical instabilities triggered at a sub-grid scale whenever the grid is 

'■^ ' coarsened. This prompted a number of studies aimed at improving the LBGK method, among which we mention the 

^ ■ unconditionally stable entromc LBGK model lg|, a family of the matrix lattice Boltzmann models j6l- [ll|, th e recent 

O , multi-step kinetic models |12l4l4| . and kinematically complete LBGK models on higher-order lattices [15[. While 

enhancing the stability of the standard LBGK model, the above approaches also have to answer questions about 

accuracy and computational efficiency. In particular, the idea of using separate relaxation times for various non- 

^-H ■ hydrodynamic modes can be realized in various ways, and some of the realizations may severely affect the accuracy of 

^ \ the simulation (see, e. g. |16l. Il7| and references therein). At present, the mainstream of lattice Boltzmann research 

^j^ ■ remains with the LBGK, due to its unsurpassed simplicity, computational efficiency and acceptable accuracy. 

!r^ Under such a state of affairs it appears that an enhancement of the standard LBGK model with respect to stability, 

»«f-\ . but without a compromise on computational simplicity and accuracy, is needed. In this work, we aim at precisely 

this kind of enhancement of the standard LBGK. Below we shall introduce a lattice Boltzmann model taking into 

account different relaxation rates for various non-hydrodynamic modes in a co-moving reference frame. This model 

assumes a LBGK form with a generalized equilibrium explicitly obtained (here, for the two-dimensional case) and 

it does not incur any significant computational overhead with respect to the standard LBGK. Extensive simulations 

of a chosen benchmark flow (roll up of a shear layer) reveal that the enhanced LBGK model features significantly 

increased stability, while retaining the accuracy of the LBGK. On the practical side, the present numerical algorithm 

is simple, requiring just a few line changes in existing standard LBGK codes. 

\—{ ' The outline of the paper is as follows: In section HIl details of the construction of the enhanced LBGK model in two 

. . 1 dimensions are presented. In section Hill the stability and accuracy of the present scheme is assessed in a benchmark 

simulation of shear layer vortical flow. Finally, some conclusions are drawn in section IIVI 
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II. ENHANCED LBGK MODEL IN TWO DIMENSIONS 



A. Moment representation 



For the sake of presentation and without any loss of generality, we consider the popular nine- velocity model, the 
so-called D2Q9 lattice. The discrete velocities are constructed as a tensor product of two one-dimensional velocity 
sets, V(^i) = i, where i = 0, ±1; thus W(ij) = (^(i) j ''^O") ) in the fixed Cartesian reference frame. Populations are labeled 
accordingly, fuj)- We start with the moment representation of the populations. To this end we recall that any 
product lattice, such as the D2Q9, is characterized by natural moments (cf. e.g. [lal)- For D2Q9, these natural 
moments are pMpq, where p = {fuj}) is the density, and 



pMpq = {f{i,j)'"li)'"u))^ P^l'^ {0,1,2}. 



(1) 



Notation (...) is used as a shorthand for summation over all the velocity indices as displayed. In the sequel, we use 
the following linear combinations to represent natural moments ([T]) 



Mqo,u^ = Mio, uy = Moi, T = M20 + A/02, N = M2Q-M02, ^xy - Mn, Q^^yy = A/12, 



Af2i, A = M22. (2) 



For the sake of completeness we recall that these are usually interpreted as the normalization to the density (Mqo = 1), 
the flow velocity components {u^^ Uy), the trace of the pressure tensor at unit density (T), the normal stress difference 
at unit density (N), and the off-diagonal component of the pressure tensor at unit density (n^-j,). The (linearly 
independent) third-order moments {Qxyy, Qyxx) and the fourth-order moment {A) lack a direct physical interpretation. 
However, as we shall see below, they are of special importance for achieving better performance of the lattice Boltzmann 
schemes. 

With the set of natural moments ([2]), populations arc uniquely represented as follows (cr, A = {— 1,1}): 



f(^ofi)=p{l-T + A). 



/(o-, 



0) 



\p{kT + N) + , 



^^xyy -^ 



/(0,A) = l^P {l^{T - N) + \Uy - XQya:x - A 

f(a,X) = TP (^ + (cr)(A)nj;y + nQ,;yy + \Qyxx) 



(3) 



We note in passing that a subset of the populations ([3]) specified by the closure relations, A = M20M02, Qxyy = 
MiqAIq2 and CX^x = MQ1M20 gives an example of a fully factorizcd population termed unidirectional quasi-equilibrium 
(UniQuE) in |18j , and it is used as an intermediate quasi-equilibrium in some constructions [1^, |l3| . We shall use a 
different route here. 

In order to construct an analog of the standard LBGK model, we further introduce central moments of the form 



pMj 



pq 



{{V(i)-UxY{V(j) -Uyff), 



(4) 



and use identity 










^xy = ^xy +UxUy, 




N^N+{ul~ul), 




T = f + u^, 




Qxyy = Qxyy + 2UyI[xy - -U^N + -U^f + UxU^, 




Qyxx = Qyxx + 2ujlxy + -^UyN + -Uyf -|- MyM^, 




A^A + 2 


^x^xyy \ ^^y^yxx 


+ AUxUyJIxy + -U^f 



(5) 



--{ul-ul)N + 



2 2 

U^Uy. 



We remark in passing that the mapping of natural moments onto central moments is nonlinear (it explicitly depends on 
the powers of the velocity components). Therefore, implementation of their relaxation in the framework of the matrix 
model with a fixed transformation matrix from moments to populations becomes involved, which may negatively 
affect efficiency and accuracy ^M,- 



Using the central moments representation, Eq. ([3|) is rewritten upon substituting ([5]) into ([3]) and rearranging 
terms, 



/(o,o) = P (l + ulul - u^) 

+ P AUxUyflxy 



9 9 



N 



u'-2 



1 -\- ZUx^xyy \ '^Uy^yxx i ^ 



.f(afi) = l^{ul+ Cru^(l - ul) - ulul) 



P 

2 

P 
2 



1 + aux + ui-Uy 



N - {2aUy + 4:UxUy)lIxy 



1 — aux — u^ 



1 (cr + 2Ux)Qxyy ~ 2UyQyxx ~ ^ I i 



2^ 



/(o,A) = ^K + Auj;(l - ul) - ulul) 



P 
2 

P 
2 



— 1 — \uy + ui ~ ul 



N — {2XUy + 4:UxUy)lix 



VJ'-'^xy 



1 — Au„ — U 



1 yA 'T' ^Uy JKc^yxx ^Ux^xyy -^ 



/(<T,A) = -{(rXuxUy + auxul + Xuyul + ulul) 



+ J I {4:UxUy + (ct)(A) + 2auy + 2Xux) Uxy + 



-ul + ul ~ aux + Xuy 



N 



u^ + aUx + Xuy 



f+{a + 2Ux)Qxyy + (A + 2Uy)Qyxx + A 



(6) 



Thus, any population on the D2Q9 lattice is uniquely represented by a linear combination of higher-order central 
moments with the coefficient of the linear combination being nonlinear functions of the flow velocity. While quite 
straightforward, this representation is helpful for the next steps of the construction. 



B. Enhanced LBGK scheme 



At the equilibrium, the higher-order (central) moments assume the following values (as dictated by the Maxwell- 
Boltzmann distribution, cf. e.g. |18l|): 



n:^, = ^°^ = Q7yy = Qllx = 0, f^"" = 2c^ A^^ = 4, 
where Cg is the speed of sound squared (reference temperature) of the D2Q9 lattice, 

1 



Let us introduce four relaxation parameters. 



'■^5^ 



W, Wfc, W3, LOA, 



and three ratios. 



ujb wa W4 

n = — , r^^ — , r4 = — . 

to UJ U! 



(7) 



(8) 



(9) 



(10) 



We consider a four-parametric family of lattice kinetic equations, written in the LBGK form (i.e. diagonal in the 
population representation). 



f{x + V,t+l)-f{x,t) = -Lu{f-r), 

where the function /* (generalized equilibrium) is constructed as follows: 



(11) 



• For any higher-order eentral moment M, introduce a hne segment connecting the current value M with the 
equihbrium value thereof, M'"^. This linear function will be parameterized with the parameter r, and denoted 

as Lr[M , M'^'i]: 

Lr[M, M""^] = (1 - r)M + rAP"^. (12) 

• In the central moment representation, Eq. ([6|), replace the second-order moments Hxy and A^ (responsible for 
the shear) by their values at the equilibrium, and replace the rest of the higher-order moments (T. responsible 
for the compressibility, and the third- and fourth-order moments, Qxyy, Qyxx, and A, respectively) by the linear 
combinations ((T^ with the parametrization according to the ratios of relaxation rates ^TU\\ . That is, if the 
shorthand notation is used for Eq. ©, 

/ = f{P, U, Uxy, N, f, Qxyy, Qyxx, A), 

we set the generalized equilibrium in the LBGK-like kinetic equation (|lip as follows: 

r = f{p,U,IlllN'^^Lrdf,f'^%LrAQ.yy,QllylLrAQy.x,Ql'i^ (13) 

In the expanded form, the generalized equilibrium (J13p reads: 

+ P ((1 - r3)[2uxQxyy + 2uyQyxx] + [(1 - r4)A + r4C^]) , 
/(Vo) = f ("' + ^"-(1 - ^D - «) + P ( ^ '"""/"" ) [(1 - n)f + 2rbc2] 

- ^ f (1 - ^3)[(cr + 2Ux)Qxyy + 2UyQyxx\ + [(1 - r^)A + T^c'S) , 

2 (14) 

/(o^A) = f ("', + A-,(l - ^^I) - ulu^^ + p ( ^~^";~"' ) [(1 - r,)f + 2r,c,2] 

- ^ ((1 - ^3) [(A + '2Uy)Qyxx + 2UxQxyy] + [(1 - ^4)^ + r4C^]) , 

t* Pf \ I 2 , , 2 , 2 2n , /u^ + CTU;r + AUj;\ .,_, N^ , o 21 

/(<T,A) = ^(crAuj;Uj^-f cru^Uy4-A%u^ + u^Uy)+pl ^ I [(1 - r6)T -h 2r(,CsJ 

+ ^ ((1 - '^3)[(ct -t- 2Ux)Qxyy + (A + 2My)Qj,^^] + [(1 - n)! + Ucj]^ . 

Here we have taken into account the actual values of the higher-order moments at the equilibrium ([7]). Using the 
standard multi-scale (Chapman-Enskog) analysis, it can be shown that ((TT|) recovers the isothermal Navier-Stokes 
equations at reference temperature To = Cg, 



dtp + daipua) = 0, 

dtUa + UpdpUa + -da{clp) dp 

p p 



Vp ( daUp + dpUa - —SapdjUj 



2 (15) 

—9a {£.pdju^) = 0, 



where _D = 2 is the spatial dimension, and where the two viscosity coefficients, the kinematic (shear) viscosity i> and 
the bulk viscosity ^ are 

Below, we shall consider the case oj = uJi, {ri, = 1) (shear and bulk viscosities equal). Moreover, Eq. p^ can be 
further simplified by neglecting all the terms of 0{u^) and higher: 

f{0fi) = P {l - 2c2 - (1 - cl)u^ + (1 - r3)[2uxQxyy + 2uyQyxx] + [(1 - ri)A + r4C^]| , 

/(Vo) = ? 1(1 - Cs)^"x + ul + (1 - u^)cl - (1 - r3)[(a + 2ux)Qxyy + 2uyQyxx] - [(1 - ri)A + r^ct]] , 

(17) 
/(*o,A) = f {(1 ~ '^s)Awy + ^2 + (1 „ ^2)^2 _ (^ _ ^^)p ^ 2uy)gy,, + 2ii,0,yy] - [(1 - ri)A + r4c4]} , 

f{a,X) = 4 { (^""^ + AUy)Cg -I- aXUxUy + U^cl + (1 - 7-3) [(cr + 2Ux)Qxyy + (A + 2Uj/)(3^y^2;] -(- [(1 - r4)A + r4cf]| . 



It can be readily checked that, by setting r^ = r4 = 1, and using the value of the speed of sound as given by Eq. 
(|S]), the generalized equilibrium ([T7|) becomes the standard second-order equilibrium of the LBGK model, whereas at 
Tb ~ r^ ~ r^ — 1, function (|T4| becomes the Maxwell equilibrium on this product- lattice (cf. [la|). 

In summary, the four-parametric enhanced LBGK model is fully explicit, and is defined by the four parameters 
(w, and the ratios rf,, r^ and r^) in the generalized equilibrium populations /* (J14p . Note that this is the maximal 
parametrization which correctly takes into account symmetries of the moments. The proposed model is readily 
implemented, similar to the standard LBGK itself. Below we shall demonstrate the gain of the present model [with 
the restricted set of parameters, Eq. (|17p ] with respect to the standard LBGK by considering a benchmark simulation 
of a shear flow. 

III. RESULTS 
A. Stability 

We shall first assess the stability of the present scheme with respect to the standard LBGK. For this purpose a 
perturbed double periodic shear layer flow is used, with initial conditions 

C/tanh(K(|- - i)) ,v < #, 



\ C/tanh(Ai(|-f 
S sin (2TT (x +\)) 



^^,f))>y>T. (18) 



as studied by Minion and Brown [l9|. L is the number of grid points in both a:; and y directions, and periodic 
boundary conditions are applied in both directions. Varying the parameter k alters the width of the shear layers, 
and this is fixed at k = 80 throughout the following. The velocity perturbation in the y-direction initiates a Kelvin- 
Helmholtz instability causing the roll up of the anti-parallel shear layers. The parameter S controls the size of the 
initial perturbation and is fixed here at 6 = 0.05. U determines the magnitude of the initial x-velocity. 

Stability regimes for the parameters r^ and r^ are considered at a fixed Reynolds number, and then the stability 
limits of the Reynolds number are considered by independently varying r^ and r^. For this the Reynolds number is 
defined as 

Re-^. (19) 

oj 2 

U is held constant throughout, with Reynolds number (at a fixed L) being varied by w alone. As U is constant, the 
initial Mach number, given by 

Ma = C/^/3, (20) 

is the same for each simulation. U ~ 0.04 is used, giving Ala w 0.07. 

To determine stability, simulations were run for a large number of time steps, T, with instability being determined 
by any deviation in the total mass inside the domain. These time steps correspond to a time given by 

TU 
t^-. (21, 

Here T = 200, 000 was used, giving t = 62.5. 

For the stability analysis a fixed grid size of 128 x 128 was used. At this grid size it was found that the standard 
LBGK scheme becomes unstable at w = 1.99692 which corresponds to Re « 20 x 10^. A series of simulations with 
the present model were then run at the same value, cj = 1.99692, with various values of the two free parameters 
rz and r^. Results are presented in Fig. [l] where points show limiting values around which many simulations with 
varying r^ and r4 were run, the points representing the limits of stable simulation. The stability domain inside these 
points corresponds to the successful (stable) simulations. The value of a; was then increased to oj = 1.999, giving 
Re « 61 X 10'^, the result also being shown in Fig. [TJ Clearly an almost convex domain of stable values is found, 
within which all combinations of ry, and r4 produce a stable result. This domain is large at the standard LBGK 
stability limit, becoming smaller as the Reynolds number is increased into values unstable for the standard LBGK. 
The accuracy of solutions within such a domain are analysed in the following section. 

Fixing a;4 = w, the stability limit of w was determined over a range of values of Ws, the results being shown in Fig. 
[51 Also shown are the results of fixing u^ — uj, and finding stable values of oj over a range of oJi. In the first case 
the highest a; to remain stable was uj — 1.99742, simulated at wa = 1.5, which corresponds to a Reynolds number of 
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FIG. 1: rs vs n stability region, on a 128 x 128 grid at the limit of stability of the standard LBGK, u = 1.99692, i?e « 20 x 10^, 
(+, solid line), and at cj = 1.999, Re « 61 x 10^ (unstable in LBKG), (x, dashed line). The region inside each line is stable up 
to at least 200 x 10^ timesteps (i = 62.5). 
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FIG. 2: Stability limit of a; on a 128 x 128 grid, for (left) varying uj^ with fixed ma = oj and (right) varying uja with fixed 
tija = uj. Values between the trend lines (solid lines) are stable. The dotted lines show the LBGK values. Significant increases 
over the standard LBGK in the stable values of uj, and therefore Reynolds number, are observed. 

about Re « 24 x 10'^. In the second case the highest a; to remain stable was uj = 1.99914, simulated at lo^ ~ 1.99125, 
which corresponds to a Reynolds number of about Re « 71 x 10"^. Varying W3 and a;4 together gives significant further 
increase in stability, for example 0-13=^4 = 1.98 remains stable up to a; = 1.999942, corresponding to a Reynolds 
number of Re w 1 x 10^, 50 times greater than the maximum stable Reynolds number using the standard LBGK. 

The grid chosen for this set of numerical experiments is too coarse to assess the accuracy of the method. It is 
well known that insufficient resolution in the present benchmark results in spurious vortices which contaminate the 
simulation. Many conventional numerical methods, as studied in Minion and Brown [19i |. are shown to produce 
spurious vortices on 128 x 128 grids at Reynold number of 0(10^). 

It should be stressed that the instability of the standard LBGK was indeed triggered by these spurious vortex 
structures, that is, due to lack of resolution. The present model is able to sustain under these circumstances even 
at much higher Reynolds numbers. We therefore proceed with the accuracy study of the present model under grid 
refinement. 



B. Accuracy 



In order to assess the accuracy of the present scheme, the Reynolds number was initially fixed at Re = 30 x 10"^, 
while the grid was doubled in each direction. The increase in the resolution stabilized the standard LBGK, however, 
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FIG. 3: Vorticity field at t 
Enhanced LBGK with ujs — 
with the Enhanced LBGK. 



= 1, on a 256 x 256 grid with Re = 30 x 10^, for the standard LBGK (left), and for the present 
uj, ixii = 1.99 (right). The spurious vortices present for the standard LBGK are completely removed 



Method 


Minimum Resolution 


LBGK 

tiJa = 1-5 

uji = 1.99 

u;3 = 1.8, aj4 = 1.95 

^3 = 1.0, a;4 = 1.90 

i03 =0.5,a;4 = 1.80 


288 X 288 
280 X 280 
248 X 248 
192 X 192 
168 X 168 
144 X 144 



TABLE L Minimum resolution for which spurious vortices are not observed at Re = 30 x 10'^ 



the spurious vortices were still present, as shown in Fig. |3l The simulation was run with the present scheme using 
the following sets of values: 1) wa = 1.0, 014 = ^, 2) ^3 ~ uj,uj4 ~ 1.99, and 3) ws = 1.0, u;4 = 1.9, as shown in Fig. 
[21 It can clearly be seen that for UJ4 ~ 1.99 the spurious vortices are completely removed. For W3 = 1.0, a;4 = 1.9 
the spurious vortices are removed even on a 128 x 128 grid. On the smaller grid the lower resolution thickens the 
shear layer slightly, however it is worth noting that on the larger grid the spurious vortices are completely removed 
for 0^4 = 1.99 without any noticeable increase in the width of the shear layers, as would occur if spurious vortices were 
suppressed through an additional artificial viscosity (cf. e. g. Ref. |19j). 

The improvement provided by the present method can be more clearly seen in Table [T] Here the approximate grid 
resolutions at which spurious vortices disappear are given for both the standard LBGK and the various setups of the 
present scheme. This again shows the present method providing a clear advantage over the standard LBGK. The 
minimum grid resolution required to remove spurious vortices on the standard LBGK is 288 x 288, compared with 
only 144 x 144 using the present method. In terms of computational time, looking at grid size alone this represents 
an eightfold reduction to give a solution of equal accuracy, although a small overhead is incurred with this method. 
This also represents a fourfold reduction in required memory. 

While the removal of spurious vortices provides obvious qualitative improvements in results, it is useful to compare 
the convergence of the current method with that of the standard LBGK. For this a solution at a high resolution, 
1440 X 1440, is produced in each case at Re = 20 x lO'^. The average differences in the x-component of velocity of this 
result compared with those of varying grid resolution are shown in Fig. |31 All setups are seen to have approximately 
the expected second-order convergence, however in each case this deteriorates below a certain resolution. In agreement 
with the results of TablclH this happens for a higher grid resolution in the standard LBGK case, and decreases to lower 
grid resolutions with the present scheme, the deterioration in second order convergence being due to the formation of 
the spurious vortices. These results also confirm that the present scheme follows the same behaviour as the standard 
LBGK, the convergence rates being equal to those of LBGK. This lends further evidence to the present scheme being 
an improvement on the standard LBGK, without compromising its underlying quality. 
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FIG. 4: Errors between solution on a 1440 x 1440 grid and varying grid sizes, showing convergence of the standard LBGK {</) 
and the present method with 1) ljs = 1.94, a;4 = cu (x), 2) lj4 = lu,lu4, — 1.97 (•), 3) ljs = 1.7,a;4 — 1.6 (+). In each case p, the 
rate of convergence, is 2.16. The solid line shows second order convergence, p — 2. (Points overlap at grid size > 240 x 240.) 
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FIG. 5: Errors in the a;-coinponent of velocity, as a percentage of average 3;-velocity, between the Enhanced LBGK on a 
240 X 240 and the standard LBGK on a 1440 x 1440 grid, for cos at fixed 014 = to (left) and 014 at fixed UJ3 = to (right). 



The same high resolution solution is used to assess the accuracy of points within an wa vs W4 domain at fixed w, 
as is observed in Fig. [1] Here a 240 x 240 grid is used as no spurious vortices are observed at this resolution at 
Re = 20 X 10'^, throughout the stable domain. Results are compared with the 1440 x 1440 LBGK solution. Plots of 
errors in x-velocity are made for three slices through this domain: 1) tUi = lo, W3 varied, 2) W3 = cj, 1x14 varied, and 
3) i03 = uji varied, with results given in Figs. [5] and |6l At this Reynolds number the average error in the standard 
LBGK case is 0.32%. It can be seen that errors are very similar to those in the standard LBGK case, throughout the 
stability domain. 

For varying ^3 alone, errors are smaller than for the standard LBGK, as is also the case for varying a;4 alone. For 
varying W3 = uj^, while for some values errors are slightly higher than for LBGK, they are of the same order. The 
increased stability of the present method has not affected the underlying accuracy of LBGK. In addition, a range of 
values of the parameters in the present model give the same quality of solution. 
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FIG. 6: Errors in the a;-component of velocity, as a percentage of average x-velocity, between the Enhanced LBGK on a 
240 X 240 and the standard LBGK on a 1440 x 1440 grid, for varying ujs = UJ4- 

IV. CONCLUSION 

An enhanced LBGK model has been developed which demonstrates a significantly enlarged domain of stability 
compared with the standard LBGK. This enhancement is obtained without any significant computational overhead 
above the standard LBGK. In that respect the present formulation can be preferable to other methods of realization of 
the additional relaxation of non-conserved modes. A large overall gain in stability is found in a benchmark simulation, 
with up to 50 times increase in accessible Reynolds number reported. This increase in stability is obtained without 
the use of artificial viscosity. A domain of relaxation parameters has been found that allow stable simulation, within 
which the accuracy of solution is consistent regardless of the choice of parameter values. This parameter range shrinks 
with increasing Reynolds number which indicates that existing approaches based on the additional relaxation times 
have a limitation as they do not provide unconditional stability. This study shows that the present approach is useful 
as it requires lower resolution grids to produce the same accuracy of solution as the standard LBGK. This enables an 
overall eight times reduction in computational effort as compared with the standard LBGK. The three-dimensional 
realization of the present model is straightforward and will be addressed in our follow-on work. 

Finally, we reiterate that the method developed herein follows the idea of arranging for different relaxation rates 
of central moments in a co-moving reference frame, as was first expressed by Geier et al [10[. However, the present 
realization is different from the cascaded LB method of Ref. [101 ■ While the cascaded LB method ]1^ stabilizes the sim- 
ulation with the help of artificial viscosity (which can sometimes result in various artifacts such as a re-laminarization 
reported in [131); the present realization retains the accuracy of the standard LBGK, as was demonstrated in the 
simulations above. In that respect, the cascaded LB method can be regarded as a sub-grid model whereas the present 
realization is more suitable for the direct numerical simulation. 
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